Method for analysing and processing seismic reflection data for the determination of a high resolution spatial velocity field for hyperbolicity correction

ABSTRACT

A method of analysing and processing seismic reflection data for the determination of a high resolution spatial hyperbolicity correction velocity field. 
     The method is characterised in that: 
     for the set of normal moveout-corrected gathers, the maximum values of the positive and negative residues of the NMO correction are determined: a time range analysis is determined located on either side of a time t 0  and of which the width is equal to not more than twice the absolute value of the maximum residual moveout, 
     a family of 2n+1 residual correction hyperbolas or parabolas are constructed, each having its apex centered on said time t 0  and, at the value of the maximum offset, presenting a value of time that is equal to one of the 2n+1 equidistant time values predetermined on the analysis range, and including the value t 0  and the extreme values of said analysis range, 
     2n+1 sets of static corrections are determined for each of the offsets, defined by the time differences presented relative to said time t 0  on the 2n+1 residual moveout correction hyperbolas or parabolas, 
     for each gather, and successively to each of the traces and for each hyperbola or parabola, the static correction associated with the offset trace is applied, and 
     the statically corrected traces are stacked together in order to obtain a set of 2n+1 stacked traces characterising correction velocity field.

BACKGROUND OF THE INVENTION

In petroleum exploration, a seismic reflection survey is a common methodfor obtaining a seismic image of the subsurface. In this method, usingappropriate energy sources, called emitters, acoustic waves aretransmitted, travel in the subsurface to be explored, and are reflectedon the different reflectors which it contains. The reflected waves arerecorded, as a function of time, on adapted receivers disposed on theground surface or in the water. Each recording (or trace) given by areceiver is then assigned to the location of the point which is situatedat the middle of the segment connecting the source to the receiver. Thisoperation is referred to as common midpoint gather.

A seismic prospecting technique, which is now conventional, is multiplecoverage. In this technique, the sources or emitters and receivers aredisposed on the ground surface in such a way that a given midpointgathers several recordings. The series of recordings associated with thesame midpoint forms what is generally called a common midpoint gather ofrecordings or traces. The set of gathers is associated with a series ofdifferent midpoints preferably located along the same line at thesurface. Based on these gathers, seismic processing serves to obtain aseismic image in the vertical plane passing through all these midpoints.The arrival time of a recorded wave varies with the angle of incidenceθ, which is the angle between the normal to the reflector at thereflection point, called the mirror point, and the direction of theincident (descending) wave. For a given gather and a given mirror point,this angle varies for each recording as a function of the offset x ofthe receiver relatively to the midpoint. Making the conventionalassumption of a homogeneous and isotropic subsurface, in plane andparallel layers, the reflections associated with each of the subsurfacereflectors, observed on a common midpoint gather, are theoreticallyaligned along hyperbolas centered on the vertical to the midpoint andcalled time/distance curves. In order to build the stack of therecordings of each gather, it is necessary to apply a correctiondepending on time, called the normal-moveout correction, which is aimedto straighten the hyperbolas to bring them theoretically to thehorizontal. Conventionally, the normal-moveout correction made is acorrection based on the following equation: ##EQU1## where: x is theoffset,

t₀ is the zero-offset source/receiver reflection travel time,

V, which is a function of time, is the average wave propagation velocityin the subsurface, and

t is the travel time after reflection associated with a source/receiverpair for offset x.

To make satisfactory NMO corrections, it is necessary to know thevelocity distribution V(t) at each midpoint. To achieve this, velocityanalyses are made at given locations on a limited number of commonmidpoint gathers. The results are then subjected to a doubleinterpolation, in time for each of the analyses, each analysis onlygiving a maximum of some twenty velocity values associated with the samevertical, and in abscissa, between the analyses, said analyses beingcommonly performed only every 40 to 50 midpoints on the average, and anaccount of the fact that the analysis is made manually, thus implyingrelatively long analysis times and relatively high processing costs.

The conventional velocity analysis consists in applying constantvelocities in succession to the common midpoint gathers, for themidpoints selected, to make the corresponding NMO corrections, and thenin stacking the corrected traces for each of the velocities used and inmanually retaining the velocities that lead to an energy peak of thestacked trace. The accuracy of the velocity field obtained by thisprocess is insufficient for a large number of more sophisticatedtreatments applied to the prestack traces (recordings), for examplemigration, inversion, measurements of effects of variations in amplitudewith offset (denoted by AVO for Amplitude Variation versus Offset),because these processes are distorted by the effects of the time andabscissa interpolations, by the inaccuracy of the velocity valuesselected, and by a signal distortion due to the NMO correction formula(1).

In an article entitled `Normal moveout revisited: Inhomogeneous mediaand curved interfaces`, published in the review Geophysics, Volume53,No.2, February 1988, pp.143-157, Eric de Bazelaire has developed anothermethod of velocity analysis, used for building improved stackedsections, called `Polystack` sections. This analysis consists inconstructing a document called BAP, which is associated with a commonmidpoint gather of recorded seismic traces.

The BAP associated with a common midpoint gather is another gather oftraces, each of which is the result of the stacking of the traces of thecommon midpoint gather after the application to each of these traces ofa static type of correction (independent of time) and different from onetrace to another. These static corrections, which have the effect ofshifting all the samples of a trace by the same time interval, arerepresentative of predetermined curvatures. Said static corrections aredefined by correction hyperbolas according to the following formula:##EQU2## where: x is the offset,

V₁ is the wave propagation velocity in the first subsurface layerexplored,

t and t₀ are the source/receiver reflection travel times for x offsetand for zero offset respectively,

t_(p), called the time of focusing depth, is the sum of the time t₀(time between the origin of the coordinates of the apex of thecorrection hyperbola) and the time between the origin of the coordinatesand the centre of the hyperbola, of which the asymptote is controlled bya unique velocity. In the aforementioned formula, the term in x² couldbe followed by at least one higher order term in x⁴ which is ignored.

Accordingly, the N traces of a BAP correspond to an investigation alongN different moveouts or curvatures.

The BAP associated with a given common midpoint gather explores thewhole field of possible hyperbolic curvatures for said midpoint, fromthe most concave at the bottom (low velocities) corresponding to lowpositive t_(p), to the most convex towards the top (imaginaryvelocities), corresponding to negative t_(p), and including the infinitevelocity (t_(p) =Y). Since the area scanned is very wide, each of theBAP has a large number of columns, generally more than 200.

Although it serves to obtain an accurate and continuous velocity field,the Polystack method nevertheless has a number of disadvantages. First,it is a costly method, because computer time is long due to the size ofthe BAP (more than 200 columns). Secondly, risks of instability existdue to automatic picking errors on the BAP. Thirdly, this method doesnot at all eliminate the undifferentiated multiples of the real events.

SUMMARY OF THE INVENTION

The present invention proposes a semi-automatic method for analysing andprocessing seismic reflection data, which helps to remedy theshortcomings of the previous methods mentioned, while preserving theiradvantages. In particular, the method of the invention helps to obtain ahigh resolution spatial velocity field, with better resolution, as wellas an improved stack section.

An aim of the present invention is a method of processingvariable-offset seismic reflection recordings or traces, in order toobtain a high resolution spatial velocity field and/or composite seismictraces. Using common midpoint gathers, this method requires subjectingthe traces of each of the gathers to a normal moveout correction, and itis characterised in that:

for the set of normal moveout-corrected gathers of traces, the maximumvalues of the positive and negative residues of the NMO correction aredetermined relative to a time t₀ corresponding to zero offset: a timerange analysis is then determined located on either side of said time t₀and of which the width is equal to at least the sum of the absolutevalues of said maximum values and at most to twice the absolute value ofthe maximum residual moveout,

a family of 2n+1 residual correction hyperbolas or parabolas areconstructed, each having its apex centred on said time t₀ and, at thevalue of the maximum offset, presenting a value of time that is equal toone of the 2n+1 equidistant time values predetermined on the analysisrange, and including the value t₀ and the extreme values of saidanalysis range,

2n+1 sets of static corrections are determined for each of the offsets,defined by the time differences presented relatively to said time t₀ onthe 2n+1 residual correction hyperbolas or parabolas,

for each gather of traces, and successively to each of the traces andfor each hyperbola or parabola, the static correction associated withthe offset trace is applied, and

the statically corrected traces are stacked together in order to obtaina set of 2n+1 stacked traces characterising a correction velocity field.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and other features and advantages of the invention will appearon a reading of the embodiment of the method of the invention and theappended drawings wherein:

FIG. 1 is a schematic representation of an array of sources andreceivers for the transmission and reception/recording of wavestravelling in a first subsurface layer to be explored,

FIG. 2 is a schematic representation of a set of synthetic seismictraces as a function of the offset of the receivers,

FIG. 3 represents a real gather of traces,

FIG. 4 represents the gather of traces of FIG. 3, after normal-moveoutcorrection,

FIG. 5 schematically represents a sample of a family of hyperbolas usedfor residual hyperbolicity correction,

FIG. 6 represents a field of investigation in an analysis range,

FIG. 7 schematically represents a BAP according to the invention,

FIG. 8 schematically and empirically represents a succession ofamplitude peaks selected on the BAP of FIG. 7 along the time axis,

FIGS. 9 and 10 schematically represent respectively the results of aconventional velocity analysis and a velocity analysis according to theinvention,

FIG. 11 is a schematic representation of a composite stacked traceobtained from the BAP of FIG. 7.

DETAILED DESCRIPTION OF EMBODIMENT OF THE INVENTION

In the device shown in FIG. 1, the source S₁ /receiver R₁ pair isassociated with seismic trace E₁ of FIG. 2 with minimum offset. Thesource S₂ /receiver R₂ pair is associated with trace E₂ with higheroffset, and so on until the source S_(Q) /receiver R_(Q) pair associatedwith trace E_(Q) obtained for maximum offset. The reflection eventsassociated with each interface, such as H, of the geological model, areshown on each of the synthetic recordings of FIG. 2 by a juxtapositionof responses ordered along time/distance curves, of which only four, C₁to C₄, are shown in FIG. 2. The time/distance curves C₁ to C₄ areportions of positive hyperbolas, of which the asymptotes are centered atthe center of coordinates. In the case of FIG. 2, the equation of thetime/distance curves characterising the travel time of the wavetransmitted at a point S and received at a point R after reflection atthe interface (or reflector) i is: ##EQU3## where: x is the offset,

t_(i) and t_(i0) are the wave arrival times respectively at the receiverwith x offset and for zero offset (receiver and source placed at thesame surface point),

and

V_(s) (stack velocity) is the average apparent velocity in the mediumtraversed by the wave above reflector i. V_(s) is expressed simply as afunction of the average propagation velocities in each of the layerstraversed by the wave, only in the case of plane and parallel layers.

The NMO corrections, which are applied before any stacking of the datacontained on the gathers associated with the same vertical, consist incorrecting recordings for the obliquity effects of the wave paths in thesubsurface, such as path S₂ PR₂ in FIG. 1, to bring them theoreticallyto vertical paths MPM. These corrections have the effect of bringing thetime/distance curves theoretically to the horizontal, a position forwhich the track stacking of the corrected traces is optimal for the dataon the time/distance curves.

Since direct recognition of the moveout pattern is not easy on realgathers, like the one in FIG. 3, the following indirect method is usedfor this recognition:

an assumption is made as to the shape of the time/distance curves of thereflections, such as I₁, I₂ and I₃,

among all the curves that comply with the assumption, the curves whichoptimise the stacking of NMO corrected traces are determined.

In conventional seismic prospecting, the only parameter that is variedin the propagation equation (3) is the value of the average apparentvelocity V_(s). A scanning of apparent velocities helps to define afamily of curves which comply with the assumption of plane and parallellayers. The conventional method for velocity analysis can be used toselect, from these curves, those that optimise the trace stack. Theresults of the conventional velocity analysis are obtained by the datumof some ten or more velocity/time pairs, the velocities V_(S) beingthose which are associated with the curves optimising the trace stack,and the corresponding times t₀. The velocity values V_(S) are thenlinearly interpolated in time in order to obtain a function V_(S) (t)associated with the midpoint analysed. A linear interpolation inabscissa x is then made between the results of the different consecutivevelocity analysis, obtaining as the final result a velocity fieldV(x,t).

In the method of the invention, normal moveout corrections are firstapplied by using either the velocity V_(S) obtained as previouslydescribed or a velocity V'_(S) obtained by a velocity field V'(x,t)known a priori. To do this, a time t₀ is assigned to the amplitude ofeach sample of each of the traces recorded with offset x, no longer itsinitial time t but said time t₀ which corresponds to it, defined by theequation: ##EQU4## where V_(S) is V(x,t) or V'(x,t).

If the calculated t₀ are not multiples of the sampling interval, aninterpolation is made between the successive calculated values.

Obviously the value of the process described above depends on theinitial assumption on the subsurface structure. The result of theapplication of conventional NMO corrections to the real example in FIG.3 reveals, in FIG. 4, the existence of residual moveout, for example onthe time/distance curves I'₂ and I'₃, which are larger as the realsubsurface geometry deviates from the initial assumptions. In this case,the results of velocity analysis are wrong, and the conventionally builtstack, made using the results of these analyses alone, is degraded.

The invention proposes in a second phase to reduce these residualmoveouts, in order to determine a high resolution spatial, continuousand better resolved velocity field, as well as an improved stack.

It is first essential to define the equations which are used tocalculate the residual corrections to be applied to the time/distancecurves.

The optimal treatment consists in applying the so-called `with delay`equation denoted P SCAN, based on the laws of geometric optics: ##EQU5##where: t₀ is the time separating the center of coordinates O from theapex A of the hyperbola,

t_(r) denotes the time separating the center of coordinates O from thecenter of symmetry of the time/distance curve (point D for hyperbola51),

V₁ is the average velocity or better, the lowest velocity of the firstmedium traversed and,

x is the offset (FIG. 5).

This equation defines a family of hyperbolas that is better adapted thanthe family of hyperbolas defined by the conventional equation (1) to thegeological media of any geometry presenting a cylindrical symmetry (orperpendicular to the profile axis).

For a given hyperbola moveout, the variation of the parameter t_(r)results in a translation of the hyperbola along the time axis, and, inseismic prospecting, corresponds to the application of a staticcorrection. By making the change of variable t_(p) =t₀ +t_(r) (wheret_(p), by analogy with optics, can be treated as a focusing depth), thepreceding equation becomes: ##EQU6## as previously indicated in (2),where 1/t_(p) -t₀ is the curvature of the time/distance curve

The normal moveout correction is supplemented according to theinvention, by a static correction which shifts by a time δt all thesamples making up each of the traces. This operation is very fast in avector computer or one equipped with a matrix processor. This makes itpossible to perform the operation, no longer as for the dynamic velocityanalyses every m kilometers, but with manual velocity picking and adouble interpolation, doing this automatically for all the samples ofall the traces and for a wide range of hyperbola moveouts. This novelmethod does not require additional trace stretching in comparison withthe NMO correction. Furthermore, it becomes possible to find velocityvariations, positive as well as negative. Since all the events, such asreal events, multiples, diffractions etc., are presented in the form ofhyperbolas, they theoretically receive equal processing, provided theyare deliberately sought. Multiples, which are manifested by moveoutsthat are quite different from those of the real surrounding reflectors,are eliminated by the processing according to the invention, contrary towhat happens in the methods of the prior art, as long as the analysisrange selected is sufficiently narrow.

To apply the second processing phase, it is first essential to definethe height in time (milliseconds) of a analysis range, which will beused on all the recordings processed. The approximate height of theanalysis range is empirically determined, for example by doubling themaximum residual moveout relative to the horizontal (corresponding to aninfinite moveout) observable on the time/distance corrected gathers. Atmaximum offset x_(max) in FIG. 4, for example, the maximum difference isobserved for time/distance curve I'₃ at maximum offset and has theapproximate value Δt=100 ms. The height of this range is advantageouslyset at 2Δt or 200 ms.

In a following step, the number 2n+1 of hyperbolas to be selected in theanalysis range is determined, and, for offset x_(max), the time intervalseparating the different hyperbolas analysed. For a correct datasampling, the signal theory defines a sufficient number 2n+1 ofhyperbolas covering the scan range 2Δt. In practical terms, a samplingis made by taking: ##EQU7## where K is a constant equal at most to √2/2and B is the passband of the seismic event processed. The value of theratio n is taken as an integer and from above. In FIG. 4, the value of Bis 50 Hz, K is 0.5, and 2n+1 is 21 for a analysis range of 200 ms.

The interval between hyperbolas for offset x_(max) is determined byp=Δt/n. In the aforementioned example, p is 10 ms. The values of t_(p)associated with 2n+1 hyperbolas serving for the static corrections andfor the creation of a document called BAP, according to the invention,are given by the equation: ##EQU8## where δ_(max) successively assumesthe 2n+1 multiple values of p covering the analysis range -np to +np,and x and V₁ are known from elsewhere.

Such a set of hyperbolas appears in FIG. 6. Each of these hyperbolas 61to 65, serving for the residual moveout correction to the previouslyNMO-corrected gather, is characterised by an index, called the ITPindex, ranging from -n to +n, in connection with the deviation δ_(max)which characterises it at offset x_(max) relative to the horizontal axisdefined by the value t_(p) =∞. By way of example, FIG. 6 showshyperbolas 61 to 65 characterised respectively by the ITP indexes -10,-5, 0, +5 and +10, hyperbola 61 presenting a maximum negative deviationof -100 ms for offset x_(max) and hyperbola 65 presenting a maximumpositive deviation of +100 ms for the same offset x_(max) relative tothe value taken for zero offset. These maximum deviations constitute thebounds of the analysis range. Such a family of hyperbolas displays ascharacteristic the fact that their asymptotes are parallel to eachother, as shown in FIG. 5, the apex A being the same.

In a following step, using the 2n+1 equations of the P SCAN hyperbolasdetermined in the previous step, the 2n+1 sets of static correctionswhich are applied to each of the traces of the gather are calculated. Aset of static corrections for a hyperbola, which is characterised by itst_(p), is composed of all the δ defined by the equation: ##EQU9## wherex assumes the different successive offset values characterising each ofthe traces. By way of example, FIG. 6 shows the value of δ for offset xjand hyperbola 61 with ITP index -10.

For each moveout characterised by its ITP index (or its correspondingt_(p)), and to the samples of each of the xi offset traces, theassociated static correction δ(xi) is applied, defined by equation (9)for the t_(p) considered.

This gives a statistically corrected gather, which is stacked to obtaina primary stacked trace characterised by its ITP index or the value ofits t_(p).

This operation is carried out for each ITP, thus giving a new gather ofstacked traces ordered according to the increasing values of the ITP.This gather is called BAP or BAP DELTASTACK. A BAP DELTASTACK is shownschematically by way of example in FIG. 7. This BAP comprises the 21primary stacked traces with ITP index 0, and +1 to +10, and -1 to -10.Such a BAP is built for each of the corrected gathers after they havebeen statically corrected as indicated above using the 2n+1 sets ofstatic corrections calculated for all the data processed. Amplitudepeaks can be observed on the primary stacked traces of each of the BAP.For example, primary stacked trace 70 with ITP index -3, in FIG. 7,displays three significant peaks 71 to 73 located approximately at 400,1350 and 1750 ms respectively.

The next step of the method of the invention consists in picking thesignificant amplitude peaks on the BAP. One way to do this consists, forexample, in first constructing the amplitude envelope for each primarystacked trace of the BAP. The BAP can then be considered as a grid ofamplitude values indexed in time and ITP indexes, and a 2D operator canbe applied to find the energy peak. The energy peaks on the BAP can befound by working with one dimension, building a representative trace ofthe maximum energy as a function of time, constructed by retaining, foreach time, the maximum of the 2n+1 amplitude values on the 2n+1 tracesof the BAP and the time concerned. The ITP index of the stacked trace ofthe BAP having the maximum value is stored in computer memory. Such atrace is shown schematically as a function of time in FIG. 8.

This trace 80 displays peaks 81 to 85 which are picked using a windowthat is slid along the time axis. Only the peaks which can be centeredon said window, such as peaks 81, 83 and 85 on trace 80, are selected.

Regardless of the picking mode selected, in the preferred embodimentdescribed, the matter is handled so that two peaks are not selected forthe same time on the BAP.

In the case of picking using a 2D operator, the previous condition issatisfied if the width of the operator is equal to the width of the BAP.In this type of unidimensional picking (FIG. 8), this condition issatisfied automatically because only one peak is selected per time.Advantageously, the width of the window or the height of the 2D operatoris equal to 2/B.

The determination of the peaks of the BAP and of the parameters (timeand t_(p)) associated with these peaks is the first result of the methodof the invention. These parameters permit the determination of a highresolution stacking velocity field at any midpoint, as well as thebuilding of a stacked trace which reconstructs the reflected eventsbetter than a conventional stack section which only uses the NMOcorrection phase of the common midpoint gathers.

It is shown that, as a first approximation, the total correctionvelocity V_(f) can be calculated from the following formula: ##EQU10##where: V_(S) is the initial NMO correction velocity,

V₁ is the wave velocity in the first layer of the medium explored,

t₀ is the zero offset time,

t_(p) is the focusing depth time.

Since all the parameters in equation (10) are known, it is possible tocalculate the value V_(f) for each of said peaks selected on each of theBAP, from the corresponding t_(p). Each value V_(f) is associated withthe time t₀ which corresponds to it. The number L of (V_(f), t₀) pairsvaries from one BAP to another. The total (V_(f), t₀) pairs serve todetermine a new velocity distribution for each midpoint.

In general, the method of the invention can be used to obtain, for eachmidpoint, a number of values 1 of velocities V_(f) that is larger thanthe number of values obtained by conventional velocity analysis,commonly ranging between 8 and 20.

By way of illustration, in FIGS. 9 and 10, a comparison is made of theresults of a conventional velocity analysis (FIG. 9) and a velocityanalysis according to the invention (FIG. 10), both obtained from thegather in FIG. 3.

To obtain a high resolution velocity field V", the knowledge of thevelocity distribution V_(f) for each midpoint helps to reduce thevelocity field calculation to a simple interpolation in time.

The result of the picking of peaks 71 to 74 on the BAP in FIG. 7 alsohelps to build a time section of which each of the traces, for eachmidpoint of the profile processed, is prepared as follows: the portions111 to 114 are selected of the primary stacked traces 70 centered on thepeaks 71 to 74, of predetermined time (shown in FIG. 7 by smallhorizontal dashes), less than or equal to the height in time of thesliding window used to pick the peaks in FIG. 8. A composite trace 110is generated (FIG. 11), associated with the midpoint by joining thetimes of portions 111 to 114 thus selected on the same axis, bypreserving the position of these portions in time, and by connectingthem by linear portions of the BAP space.

Alternatively, a portion of the processing consists in successivelyapplying to each trace of each NMO-corrected gather, 2n+1 sets of staticcorrections defined, for the offset corresponding to the trace, by thetime deviations presented relatively to t₀ on the 2n+1 static correctionhyperbolas. The corrected traces associated with each midpoint are thengrouped into 2n+1 groups, each consisting of the traces which haveundergone the static corrections corresponding to one of the correctionhyperbolas, and in stacking the corrected traces of each group to obtaina set of 2n+1 primary stacked traces corresponding to the 2n+1 staticcorrection hyperbolas, said set constituting the set of seismic tracescharacterising the fine NMO correction velocity field for the midpointconsidered.

In the foregoing discussion, reference has been made to families ofcorrection hyperbolas answering to formula (2). However, it is possibleto use families of parabolas instead of said families of hyperbolas. Theresolution of the equation defined by formula (2) leads to the followingformula: ##EQU11##

A series development of formula (2') leads to a new formula (2") of thetype: ##EQU12##

If x is small relative to V₁, the x⁴ and higher order terms can beignored. In this case, formula (2") is reduced to: ##EQU13## which isthe equation of a parabola.

The steps of the method of the invention can be carried out with thefamily of parabolas defined by formula (2'"). Thus, in FIG. 6,hyperbolas 61 to 65 can therefore be replaced by the correspondingresidual moveout correction parabolas on the NMO-corrected gather, andthen by proceeding as indicated in the description of the preferredembodiment of the invention.

We claim:
 1. A method of processing variable-offset seismic reflection traces to obtain a set of traces that can be used for the determination of at least one of a high resolution spatial velocity field and composite seismic traces, wherein the method comprises the steps of:(a) obtaining normal moveout (NMO) corrected seismic traces of common midpoint gathers; (b) subjecting the traces to a dynamic obliquity correction comprising the steps of: (c) determining the maximum values of the positive and negative residues of the NMO corrected gathers relative to a time t_(o) corresponding to zero offset; (d) determining an analysis range, the time-width of which is between the sum of the absolute values of the maximum values of the positive and negative residues and twice the absolute value of the maximum residual moveout; (e) constructing a family of 2n+1 residual correction curves, taken from the group consisting of hyperbolas and parabolas, having their apexes centered at the time t_(o) ; (f) determining a time interval between the 2n+1 equidistant residual correction curves using the maximum offset and the analysis range; (g) determining 2n+1 sets of static corrections for each of the offsets using the time interval relative to the time t_(o) and the 2n+1 residual correction curves; (h) applying the 2n+1 sets of static corrections to each trace of the gather in succession and in accordance with the appropriate offset; and (i) stacking the statically corrected traces and obtaining a set of 2n+1 stacked traces representing a correction velocity field.
 2. The method of claim 1, wherein the step of obtaining the normal moveout corrected traces of the common midpoint gathers comprises using conventional velocity analysis at a plurality of predetermined midpoints.
 3. The method of claim 1, wherein the step of obtaining the normal moveout corrected traces of the common midpoint gathers comprises using a known velocity field.
 4. The method of claim 1, wherein the steps of constructing the family of 2n+1 residual correction curves and determining the time interval comprises the step of separating the curves at the predetermined value of the maximum offset, wherein n is the nearest integer of the ratio of one-half of the analysis range (Δt) and K/B, and wherein K is a coefficient at most equal to (2)^(1/2) /2 and B is the passband of the common midpoint recordings.
 5. The method of claim 1, wherein the family of static correction hyperbolas are defined by the equation: ##EQU14## wherein x is the offset, t and t_(o) denote the source/receiver travel times for offset x and for zero offset respectively, (t-t_(o)) is the value of the static correction as function of the offset x, t_(p) is the sum of the time t_(o) and the time separating the origin of the coordinates from the center of the hyperbola, and V₁ denotes the propagation velocity of the waves transmitted in the first subsurface layer explored.
 6. The method of claim 1, wherein, for low values of x relative to V₁, the family of static correction parabolas is defined by the general equation:

    (t-t.sub.o)=(1/2)px.sup.2

wherein x is the offset, t and t_(o) are the source/receiver reflection travel times for offset x and for zero offset respectively, (t-t_(o)) is the value of the static correction as function of the offset x, p=1/(t_(p) V₁ ²), t_(p) is the sum of the time t_(o) and the time separating the origin of the coordinates from the center of the parabola, and V₁ is the propagation velocity of the waves transmitted in the first subsurface layer explored.
 7. The method of claim 1, further comprising the steps of:(j) selecting the peaks on each of the sets of the 2n+1 stacked traces using a sliding window with a predetermined time height, the time height being equal to about 2/B, and selecting only the centered peaks in the window; (k) associating each peak with a value of t_(o) and t_(p) ; and (1) calculating the velocity V_(f) corresponding to the time pair (t_(o), t_(p)) using the equation: ##EQU15## where V_(s) is the velocity used in the normal moveout correction.
 8. The method of claim 1, further comprising the steps of:(j) selecting the peaks on each of the sets of the 2n+1 stacked traces using a sliding window with a predetermined time height, the time height being equal to about 2/B, and selecting only the centered peaks in the window; (k) selecting portions of the stacked traces centered on the peaks which are of predetermined time-length, the time length being equal to about the time-length of the sliding window; (l) generating a composite stacked trace associated with the midpoint corresponding to the set of the 2n+1 stacked traces by placing the portions of the 2n+1 stacked traces on the same time axis; and (m) fixing the position of the portions on the time axis. 